Systematic generation of finite-range atomic basis sets for finear-scafing calculations 
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Basis sets of atomic orbitals are very efficient for density functional calculations but lack a sys- 
tematic variational convergence. We present a variational method to optimize numerical atomic 
orbitals using a single parameter to control their range. The efficiency of the basis generation 
scheme is tested and compared with other schemes for multiple C, basis sets. The scheme shows 
to be comparable in quality to other widely used schemes albeit offering better performance for 
linear-scaling computations. 
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The last few years have seen the development of jim- 
plementations of the density functional theory (DFT)El in 
which the computer time and memory scale linearly with 
the number N of atoms in the system studiedHQ. These 
so-called order- (0(A^)) methods have increased consid- 
erably the need of accurate and efficient basis sets of finite 
range. While high accuracy can be achieved with flexi- 
ble linear combinations of atomic orbitals (LCAO), high 
efficiency requires the orbitals to be as localized as pos- 
sible. Numerical atomic orbitals (NAO's) are well suited 
to linear scaling methods because they are very flexible, 
can be strictly localized, and few of them are needed for 
accurate results. Their main drawback is the lack of a 
systematic procedure to ensure a rapid variational con- 
vergence with respect to the number of basis orbitals and 
to the range and shape of each orbital. 

In the context of the ab initio pseudopotential method 
for solids, an early proposal were the 'fireballs' of Sankey 
and Niklewski: solutions of the radial Schrodinger equa- 
tion for an isolated pseudo-atom confined in a spher- 
ical hard potential boxcl. Subsequent works proposed 
different._u£cipes to find multiple-^ and polarization 
orbitalsauLl. In a recent workQ, a method was proposed 
to optimize the shape of the orbitals by substi1j*i±i»c 
the hard box by a soft confining spherical potentiaH'BlllJ. 
This confining potential, which may be different for each 
atomic orbital, depends on a series of parameters which 
determine the orbital's shape. The parameters are then 
adjusted to minimize the energy of a prototype molecule 
or solid. If the confining potentials diverge at given cut- 
off radii, the orbitals become strictly zero beyond those 
radii. However, if the cutoff radii themselves are included 
as variational parameters, without constraints to impose 
a small range, the resulting orbitals tend to become very 
extended, with long tails that generally have no partic- 
ular significance for the condensed system, but which 
limit severely their efficiency. In the present work, we 
propose a simple procedure to compress the orbital radii 



by introducing a fictitious pressure. This allows to bal- 
ance efficiency versus accuracy in a continuous and well 
controlled way. In addition, we evaluate the variational 
completeness of the resulting orbital shapes, by adding 
additional degrees of freedom, and by exploring alterna- 
tive generation procedures and comparing their relative 
merits. 

Our basis orbitals are products of spherical harmon- 
ics times numerical radial functions centered on atoms. 
The quantum chemistry literature typically distinguishes 
between core, valence, polarization, and diffuse basis or- 
bitals. In our case, core states are eliminated by the use 
of norm-conserving pseudopotentialstil. The explicit de- 
scription of semicore electrons as valence is performed 
with the same methods described here, but using a pseu- 
dopotential for which the semicore electrons occupy the 
ground state and the valence electrons occupy the first 
excited state (with a radial node). In previous works we 
have designed a specific numerical method for polariza- 
tion orbitalsQ, but here we will use the same methods for 
valence and polarization orbitals. We will not consider 
diffuse orbitals in this work. 

When several basis orbitals with the same center 
and angular momentum are used to expand the valence 
states, we follow the standard quantum chemical termi- 
nology and call them first-C orbital, second-^ orbital, etc, 
even though there are no C, exponent coefficients in our 
orbitals. We use a different method to generate the first-^ 
orbitals than that for the subsequent-C orbitals. For the 
first-C orbitals we solve the radial Schrodinger equation 
for a potential given by the sum of the full (screened) 
nonlocal pseudopotential corresponding to the angular 
momentum of the orbital, and a confining potential of the 
form V{r) ~ Vo cxp [—{vc — ri)/{r — r^)] /(rc — r) which 
depends on three parameters r^, Vq, and the cutoff ra- 
dius rc- These parameters are different for each basis 
orbital and define its range as well as its shape by allow- 
ing a depression of the tail. Other confinement schemes 



have beenproposedQclil3 and are compared with this one 
in Ref. Q. To generate the second and subsequent-^ 
orbitals we will use and compare two possible methods. 
The first one is based on the concept of chemical hardness 
(CH) and defines the difFerent-C orbitals as the deriva- 
tives of the ground-state wavefunction of the potential 
(pseudo plus confining) with respect to the charge of the 
atomO. In this scheme, there are no independent param- 
eters to fix the shape of the higher-than-first-C orbitals. 

The second scheme to generate higher-^ orbitals was 
inspired by the "split valence" (SV) method which is 
standard in quantum chemistry, where othitals are given 
by fixed linear combination of gaussiansllj. The second- 
C (or triple etc) orbitals are obtained by "splitting" the 
slowest-decaying gaussian(s) to act as independent ba- 
sis orbital(s). The SV was adapted to numerical atomic 
orbitals by constructing a double-^ orbital as one that 
reproduces the tail of the first-C from lajgatching radius 
outwards, and runs smoothly inwardsErQa. Higher-C or- 
bitals are obtained repeating the procedure at different 
radii. 

A variational optimization of a basis set constructed 
as described above can give orbitals with too long cutoff 
radii re- In order to reduce their range in a system- 
atic way we introduce a parameter P with dimensions 
of pressure (that we will call "pressure" henceforth) and 
minimize the "enthalpy" E + PV, where E is the total 
energy of some reference system and V — (47r/3) r^^ 
is the sum of the volumes of the basis orbitals . The 
convergence of calculated properties with respect to or- 
bital range is thus controlled by a single parameter, much 
in the same way as the planewave cutoff controls the con- 
vergence of a plane wave basis set. 

The reference system for which E + PV is minimized 
is a molecule or solid in which the atoms considered have 
a prominent role, and which is small enough to allow 
many selfconsistent calculations with different basis pa- 
rameters. The derivatives of E with respect to those 
parameters are generally|-not available, and we use the 
downhill-simplex methodtJ to minimize it. The basis or- 
bitals depend on the described parameters in a non-linear 
way, and several local minima are found in many cases. 
This is to be expected because different combinations of 
parameters can produce approximately the same optimal 
shape. Since our parameters have no special physical sig- 
nificance, any low local minimum is in principle equally 
acceptable, even though the multiple minima produce a 
somewhat unpleasant "noise" in the results reported be- 
low. 

Fig. shows the cutoff radii of the first-C orbitals of 
Si, Au, and Pb as a function of the pressure parame- 
ter P. The basis optimizations were performed in their 
corresponding bulk solids, with a so-called double-^ po- 
larized (DZP) basis set: in Si there are double-^ s and 
p shells and single-^ d orbitals; in Au there are double-^ 
s and d, and single-^ p; in Pb the 5d semicore electrons 
are included in the valence as double-C, as are the s and 
p shells, while the 6d have a single-C. The second-C or- 
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FIG. 1: Cutoff radii of the first-(^ basis orbitals of Si, Au and 
Pb, as a function of the fictitious pressure parameter P 



bitals were generated with the SV scheme. Polarization 
orbitals are obtained in the same manner as the other 
atomic orbitals. 

To give an idea of how the orbital radii affect the ba- 
sis efficiency. Fig. || shows the CPU time required for 
the calculation of a selfconsistent step of bulk silicon, as 
a function of the pressure P used to generate the basis. 
The accuracy of the results, as the orbitals contract, is 
addressed in Table |, which shows the variation in lattice 
parameter, bulk modulus, and cohesive energy witloj-E. 
The results were obtained using the Siesta methodQEJ, 
with a well converged real-space integration grid. They 
are compared to experiment and to well-converged plane 
wave calculations, performed with a specific pro gran i de- 
signed to use exactly the sanm,pseudopotentia]lliHla, ep. 
change correlation functionally, and fc-grid samplingEZI 
used in Siesta. The cohesive energy is calculated as 
the difference between the bulk total energy per atom 
(with the chosen basis set) and an atomic calculation in 
which the radial Schrodinger equation is solved numeri- 
cally, without any constraint to the shape or range of the 
orbitals. With this definition the cohesive energy car- 
ries the variational character of the total energy (higher 
binding energies for better basis sets). 

It can be seen that a moderate pressure of ~ 0.2 GPa 
produces a drastic reduction of the orbital radii, with 
a correspondingly large reduction of CPU time, without 
a significative change in the results. Larger pressures 
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FIG. 2: CPU time for a selfconsistency step of bulk Si (16 
atoms per cell) versus the fictitious pressure P used to com- 
press the cutofi^ radii of the DZP basis orbitals. 



TABLE I: Comparison of structural properties of different 
systems as a function of the pressure parameter P (in GPa) 
used to generate their basis sets. Lattice parameters a in A, 
bulk moduli B in GPa and cohesive energies Ec in eV. The 
bulk moduli were obtained Iw-ifitting the total energy with a 
Murnaghan equation of statalEI. A double-(^ plus polarization 
basis was used in all cases. In Pb semicore states where also 
used. 
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produce additional, though more moderate gains in basis 
efficiency, but at the expense of nonneglegible changes in 
the results. That small pressure of 0.2 GPa seems to be a 
threshold up to which only the very low, not significant, 
tails are removed. 

The relative merits of the SV and CH methods to gen- 
erate the second-C orbitals are considered in Fig. ^. For 
the SV case, two curves are plotted. In one of them, 
the inner matching radius of the sencond-C orbitals is 
optimized for every value of P. In the other ojic, it is de- 
termined by a standard automatic criterionQ, by which 
the norm of the first-^ orbital beyond the matching ra- 
dious has to be equal to a given "split-norm" parameter 
value of 0.15. Fig. ^ shows the optimized value of this 
parameter, which does not differ much from the standard 
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FIG. 3: Equilibrium lattice constant (a), bulk modulus (B) 
and cohesive energy (Ec) of bulk Silicon as a function of the 
fictitious pressure parameter P. A double-(^ plus polarization 
basis was used. The second-(^ orbitals were generated using 
the chemical-hardness (CH) and split-valence (SV) schemes. 
For the latter, results are shown for orbitals whose inner 
matching radii were generated with a constant split-norm pa- 
rameter of 0.15 or optimized variationally for each value of P 
(what resulted in the split-norm parameters shown in Fig. ^ . 



value. As a consequence, it is not surprising that Fig. ^ 
shows a similar quality of the results using the optimized 
and standard values. The quality is also similar for the 
CH method, which does not depend on any variational 
parameter. Again, this is not surprising, in view of the 
similarity of the resulting shapes of the second-C orbitals, 
which are compared in Fig. ^ to our SV orbitals and to 
a typiGal quantum-chemistry gaussian-based polarization 
orbitalt3. We may then conclude that the different gen- 
erating schemes of second-^ orbitals compared here yield 
basis sets of similar quality. Our SV scheme, however, of- 
fers higher efficiency for linear-scaling computations since 
the range of the higher-^ orbitals may be restricted to 
their inner matching rjadius, without any reduction of 
the variational freedomQ. 

Finally, we explore to what extent the orbital shapes 
generated with the described schemes differ from opti- 
mal. To this end, we have added spherical Bessel func- 
tions to our generated orbitals, not as additional basis 
functions but to change the shape of the orbitals in a 
DZP basis, introducing the coefficients of the linear com- 
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FIG. 4: Optimal value of the split-norm parameter, which 
determines the inner matching radius of the second-(^ orbitals 
of silicon generated with the split- valence scheme. 
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FIG. 5: Radial shape of the first and second p-orbitals of 
Si. The second-(^ orbital was generated using the chemical- 
hardness (CH) and split- valence (SV) approaches described in 
the text, bi addition, we show the second-^ gaussian orbital of 
Huzinagatj. The second-^ orbitals have been orthogonalized 
to the first-^ one to facilitate the comparison. 



bination as the parameters to be optimized. Table ^ 
shows the effect in the total energy for bulk sihcon. as 
subsequent Bessel functions are added to optimize differ- 
ent orbitals. The energy reduction is quite moderate, and 
considerably smaller than that obtained by introducing 



additional basis orbitals. This is true even in the case 
of the higher-C orbitals, whose shape depends on just 
one parameter. It can be thus concluded that the radial 
shapes of the basis orbitals are indeed well optimized by 
the variational freedom contained in the confining po- 
tential, and by the physically motivated schemes used to 
generate the higher-^ orbitals. 

In conclusion, we have developed a systematic method 
to construct accurate and efficient atomic basis orbitals 
for linear-scaling DFT calculations. The range of the 
basis sets is controlled by a single parameter, that al- 
lows to monitor their convergence with range in a simple 
and systematic way. By comparing different generation 
schemes, and by studying the effect of additional varia- 
tional free dom, we have found that our method produces 



Basis Size 



AE (meV) 



DZP not optimized 230 
DZP optimized 40 
DZP 4 Bessels in first ( 33 
DZP 4 Bessels in second ( 33 
DZP + F 22 
DZ2P + F 16 

TABLE II: Test of the quality of the DZP optimized ba- 
sis set of silicon. Second-i^ orbitals were generated with the 
split- valence method. The energies AE are per atom and rel- 
ative to the converged plane wave result. The F stands for 
the addition of a / angular momentum shell. The 2 in the 
DZ2P denotes the addition of a second to the d polarization 
orbital. The non-optimized basis was obtained with a hard 
potential!] (the radii are as long as in the DZP optimized case) 
and a standard split-norm parameter of 0.15. A zero pressure 
parameter P was used in all the cases. 



nearly optimal shapes in multiple-^ polarized basis sets. 
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